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We discuss algorithms which allow to calculate identical two-particle correlations 
from numerical simulations of relativistic heavy ion collisions. A toy model is used 
to illustrate their properties. 


1 Introduction 

The current relativistic heavy ion program at CERN and BNL aims at inves¬ 
tigating the equilibration properties of hadronic matter at extreme tempera¬ 
tures and densities where quarks and gluons are expected to be the physically 
relevant degrees of freedom for particle production processes. The theoret¬ 
ical discussion of these collision systems is complicated by their mesoscopic 
character. They are not sufficiently small to allow for an analytical descrip¬ 
tion in terms of elementary processes. They are not sufficiently large to take 
a description in terms of macroscopic observables for granted. Even if simple 
thermodynamically or hydrodynamically inspired models account for the data, 
the task remains to understand the microscopic origin of their success, and to 
establish to what extent such an agreement is necessary or accidental. 

Numerical simulations of relativistic heavy ion collisions are an adequate 
tool to attack these questions. They allow to propagate sufficiently many 
degrees of freedom from some initial condition, and they thus test the conse¬ 
quences of different assumptions about the microscopicdn medium dynamics. 
Many event generators have been developed to this aimEl. Irrespective of their 
particular model assumptions, they return for each simulated event a set of 
N m phase space emission points { (r», p*, together with informa¬ 

tion, e.g., about the identity (PID) of the particle i, the parent resonance it 
originates from, the number of rescatterings it has undergone, etc. 

The major aim of event generator studies is to reduce step by step the 
class(es) of dynamical scenarios consistent with experimental data, thus ob¬ 
taining an increasingly better picture of the physics. To this end, one mainly 
compares the simulated particle ratios and one-particle momentum spectra to 
experiment. However, one can imagine different dynamical scenarios which all 
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account for the same simulated momentum distributions {(Pi)}ie[i,jv m l but dif¬ 
fer with respect to the simulated space-time structure {(fj, • Iden¬ 

tical two-particle correlations C^q, K), here written as a function of the mea¬ 
sured pair momentum K = ^(pi +p 2 ) and relative pair momentum q = pi — p 2 
are sensitive to the latter and can thus play an important complementary role 
in specifying the parameter space of event generators. Here, we discuss algo¬ 
rithms of the type 


{ f =>C(q,K). (1) 

V ) ra£[l,iVev] 

For notational convenience, we restrict the following discussion to one particle 
species only, like-sign pions say. Also, we neglect final state Coulomb and 
strong interactionscl. The set {(fj, pj, U)} ie [\,N m \ denotes then the phase space 
emission points of the N m like-sign pions generated in the m-th simulated 
event. The event generator simulates thus a classical phase space distribution 

JV ev N m 

Pclass)**, P: I) = EE <5 (3) (r - hi) <5 (3) (p - pi) S(t - U) . (2) 

eV — 1 A —1 


2 From phase space densities to momentum correlations 

The standard particle interferometric analysis of heavy ion collisions starts 
from the relation between the measured momentum correlations of identical 
particles and the quantum mechanical Wigner nhase space density S(x,I\) of 
particle emitting sources in the collision region El, 


C(q,K)=Af 



1/ d 4 x S{x, K) e 1 ' 


Jd 4 xS(x,p 1 ) f d 4 y S (y, P2) 


( 3 ) 


This expression can be generalized to include final state interactions □. De¬ 
pending on the definition oLthe two-particle correlator in terms of single- and 
two-particle cross sectionsCm (respectively depending on its construction from 
experimental data), the normalization is either Af = 1 or Af < 1. The calcula¬ 
tion oLdeviations of Af from 1 typically involves multiparticle symmetrization 
effectsBa, which we do not review here. Our discussion will thus be for J\T = 1. 

In the light of Eq. calculating the two-particle correlator C(q, K) 
from an event generator output amounts to specifying the emission function 
S(x,K) in terms of the generated phase space points {( fj, p ^, ti )|zF fi.Af^i, i.e., it 
amounts to an interpretation of (r i , p,, U ). In sections |2.l| and |2.2|, we focus on 
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two different interpretations referred to as quantum mechanical and classical , 
though they are both conceptually on an equal footing. Before doing that, 
we list here three requirements which an algorithm implementing ( 0 ) should 
fulfill: 

1. Event generators use probabilities to describe the collision process. Bose- 
Einstein correlations are obtained from squaring production amplitudes. 
Hence, the generated momenta p;, p j do not show the enhancement at 
low relative momenta characteristic for Bose-Einstein final state sym- 
metrization. The prescriptions (JTJ) should calculate this effect a posteri¬ 
ori. 

2. The strength of Bose-Einstein correlations depends on the distance of the 
identical particles in phase space, not in momentum space. We thus re¬ 
quire the prescriptions (Q) to use the entire phase space information, and 
not only the generated momentum information. ‘Weighting’ or ‘shifting’ 
prescriptions which are based only on the latter may successfully match 
the measured momentum correlations but obviously do not allow to test 
the simulated space-time structure. 

3. In general, Bose-Einstein statistics can affect particle multiplicity distri¬ 
butions during the particle production process. However, event genera¬ 
tors are tuned to reproduce the average particle multiplicities in agree¬ 
ment with the data. We require hence that the prescriptions ( 0 ) conserve 
particle multiplicities. If the correlator is interpreted as a factor, relat¬ 
ing measured two-particle differential cross sections to those simulated 
by the event generator, d 6 a mea , s /d 3 p 1 d 3 p 2 = C(q, K) d 6 a sim /d 3 p 1 d 3 p 2 , 
then this implies tllH that Af < 1 in ©>■ This however plays no role 
in what follows since the space-time analysis of correlation data can be 
based entirely on the momentum dependence of C(q, K), disregarding 
its absolute normalization. 


2.1 “Quantum” interpretation of event generator output 

In the quantum mechanical interpretationBi&0, the set of phase space points 
{{(fj,pi,ti)} ie [ li jv m ]}me[i,Ar e Y] * s jWfflciated to the centers of Gaussian one- 
particle wavepackets according to Lrcffl 


(Pi,r i,U) —> /*(x,ti) = 


( 7 T < 7 2 ) 3 / 4 


-oMx-f,) 


*P * 


( 4 ) 


These wavepackets fi are quantum mechanically best localized states, i.e., they 
saturate the Heisenberg uncertainty relation with Axi = a/y/2 and A pi = 
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1 /y/2a. 

Taking only two-particle symmetrized contributions into account, i.e., adopt¬ 
ing the so-called pair approximation, all spectra can be written in terms of the 
functions s,(p), which coincide with the corresponding one-particle spectrum 
for a source emitting only one particln^t the phase space point i. For the one- 
and two-particle correlator one finds clQ 


E, 


d 3 N 1 ^ 


C(q,K) = l + e~ a 

Si(p) =7 t~ 312 a 3 e~ a, ( P ~ pi 


y^Nev 

Z—/m=l 



y^-Nev 

Z—/m=l 

^m(Pa) Vm(Pb) ~ Eill 1 s i(Pa) S*( Pb) 



(5) 

( 6 ) 
(7) 


Here, the terms subtracted in the numerator and denominator of C(q, K) 
are finite multiplicity corrections which become negligible for large particle 
multiplicities!! This algorithm is consistent with an emission function S(x, K) 


S(x, K) = J dii dpi dti PdassCf*, Pi, U) s 0 (x - f u t - K — p*), (8) 

s 0 (x,K) = ^5(t)e-^ x - a K , (9) 

which is given by a folding relation between a classical distribution /? c iass of 
wavepacket centers and the elementary source Wigner function sq ( x , K). The 
latter has at emission an optimal quantum mechanical phase space localization 
around phase space points (r pi,U) with spatial uncertainty cr and momentum 
uncertainty 1/cr. 

In the algorithm (§-©> the spectra are discrete functions of the input 
(fj, p i,ti) but they are continuous in the observable momenta pi, P 2 and hence, 
no binning is necessary. Also, the number of terms to be calculated in ®) and 
(S) grows linear in N, not quadratic as in other prescriptions. However, both 
the one- and two-particle spectra depend in this algorithm on the wavepacket 
width cr. The role of this parameter will be discussed in the context of our toy 
model study in section |j. 


2.2 “Classical” interpretation of event generator output 

In the classical interpretation!0, the phase space points (r,, p,, f) are seen as 
a discrete approximation of the on-shell Wigner phase space density S(x,p), 

S(x,p) = p c iass(x, p,t). (10) 
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This defines both the two-particle correlator via (^|) and the one-particle spec¬ 
trum via the ^-integration over ©• To simplify a numerical implementation 
of this calculation, it is useful to replace the <5 functions of p c iass in momentum 
space by rectangular ‘bin functions’!!], or by properly normalized Gaussianslm 
of width e, 

6 pL = (^2)3/2 ex P HP* - P)Ve 2 ) • (11) 

In the limit e —> 0, these Gaussian bin functions reduce to the properly normal¬ 
ized 5-functions 5^ (p,; — p) and we regain Eq. §)• The one-particle spectrum 
and two-particle correlator read then BE] 


Er, 


d 3 N 

<Pp 


N ev N rr 


d 4 a 


-i "'ev - 1 > m 

s ^p) = at EEC p , 


m =1 i= 1 


( 12 ) 


C(q,K) = l 


y^-ATev 
Z-^m =1 


y N ™ 5V 

Z—/i—1 p 


(<0 


P,,K' 


D ig u t i -iq-f i 


V ' V,, 


(Ck) : 


y^^Vev 

Z-/m=l 


(eS «, 


(0 


=1 Pi,Pi 


) (e£s «f‘> 


Pj ,P2 


- Et m !4 


5 ( - e) 

Pi Pi,P2 


The correlator © is the discretized version of the Fourier integrals in (^). 
The subtracted terms in the numerator and denominator remove discretization 
errors which would amount to pairs constructed from the same particles. The 
‘classical’ algorithm has the same important advantages as the ‘quantum’ one, 
it involves only 0(N m ) rather than 0(N^ l ) numerical manipulations per event, 
and its observables are discrete functions of the input, but continuous functions 
of the measured momenta; hence, no binning is necessary. ‘Classical’ and 
‘quantum’ algorithm then differ only with respect to two points: 


1. The classical algorithm has no analogue for the Gaussian prefactor exp ( 
— a 2 q 2 /2) in (|6|) which is a genuine quantum effect stemming from the 
quantum mechanical localization properties of (^). 

2. For the choice a = 1/e, the bin functions 5E p is the classical counterpart 
of the Gaussian single-particle distributions Sj(p). Finite event statistics 
puts a lower practical limit on e, but the physical momentum spectra 
are recovered in the limit e —> 0, i.e., e should be chosen as small as 
technically possible. In contrast, the quantum algorithm treats a as 
finite physical localization of the single particle wavepackets. The limit 
a —> oo (which corresponds to e —> 0) is not the physically relevant case. 
It amounts according to (j|) to an emission function which is infinitely 
extended in space, and hence 0 lirn^oo C(q, K) = 1 + 5 q .o- 
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3 The Zajc model 


Our final aim is to apply the above algorithms to event generators of relativistic 
heavy ion collisions and high piergy elementary processes. A first step in this 
direction will be reported in 113 . Here, we illustrate the properties of these 
algorithms for a toy model. 

The strategy is to start from a simple expression for the classical distri¬ 


bution pciassj from which the results of the algorithms of sections 2.1 and 2.2 


can be calculated analytically. In a second step, we generate then in a Monte 
Carlo procedure sets of phase space points {(fj, p*, i»)}*e[i,iv m ] according to 
this distribution, and we test the numerical results against our exact analyt¬ 
ical expressions. This allows us to make in an analytically controlled setting 
statements about i) the cr-dependence of the quantum algorithm, ii) the e- 
dependence of the classical algorithm (especially: how small e has to be chosen 
to extract the limit e —> 0 numerically) and iii) the statistical requirements for 
the algorithms to work. 

Here, we study these questions for the Zajc model which describes.^, 
position-momentum correlated phase space distribution of emission points E3 


Pci a a J ss( r ’ p> t) = K exp \ - 


Af s — E„ 


1 


2(1-s 2 ) 
N 


_ _ r ' P 

Rq RqPq 


p 2 

M) 


(2t r) 3 P 3 P 3 (1 - s 2 ) 3 / 2 ' 


m, (i4) 

(15) 


This distribution is normalized such that the total event multiplicity calculated 
from the one-particle spectrum @ is N. The parameter s smoothly inter¬ 
polates between completely uncorrelated and completely position-momentum 
correlated sources. For s —> 0, the position-momentum correlation vanishes 
and we are left with the product of two Gaussians in position and momentum 
space. In the opposite limit 



(16) 


the position-momentum correlation is perfect. The corresponding phase space 
volume of the distribution is proportional to 


V = 2 3 Rq P 3 (1 - s 2 ) 3/2 , 


(17) 


and goes to zero for s —> 1. This strong s-dependence allows to study the 
performance of our numerical algorithms for different phase space volumes. 
In the following subsections, we discuss the s-dependence of the one- part icle 
spectrum and two-particle correlator, focussing in sections 3.1 and 3.2 on 


analytical results, and comparing these in section 3.3 to numerical calculations. 
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3.1 The Zajc model in the classical algorithm 


Starting from an analytical emission function S(x, p) = p c iass(x, p, t) according 
to the one-particle spectrum and two-particle correlator read 


dN N ( p 2 

Ep d ^ p (2 t r) 3 / 2 P 0 3 6XP V ' 2P* 

C(q, K) = 1+ exp ( - q 2 P 2 lass C 

^class( e ) = -^0 (1 — s 


1 + c 2 /2Pq (2P 0 R 0 ) 2 1 


2 /2-Po 


(18) 

(19) 

( 20 ) 


We recover the physical HBT radius from (^) in the limit e —> 0 or by inserting 
© directly into !)• The remarkable fact is that for sufficiently large slid, 


S > Scrit " ^ (2i? 0 Po) 2 ’ (21) 

the HBT radius P 2 lass becomes negative and the two-particle correlator shows 
an unphysical rise of the correlation function with increasing q 2 . The change of 
sign in (|l9|) seems to be related to the violation of the uncertainty relation by 
the emission function ©: In the limiting case Z 2 class reduces to a perfect 

position-momentum correlated source which is not consistent with the Heisen¬ 
berg uncertainty relation. Also, the phase space volume V in (|l7j) drops below 
1 for s > s cr it. In this sense, the distribution Pciass I s a quantum mechanically 
allowed emission function 5(a;, p) = p c i aS s(x, p,t) for s < s cr ;t only. 

Clearly, an ansatz for a quantum mechanical Wigner phase space density 
S(x,I\) has to satisfy quantum mechanical localization properties. However, 
realistic event generators satisfy this requirement anyway. The practical im¬ 
portance of (this limiting case of) the Zajc model is thus that it provides an 
(extreme) testing ground for our numerical algorithms. Analytically, we con¬ 
clude already from (^o|) that a bin width sufficient to investigate the limit 
e —» 0, has to be small on the scale of the width of the generated momentum 
distribution, 

e<^V2P 0 . ( 22 ) 

At least for the present model study, this requirement is independent of the 
strength of position-momentum correlations in the source. 
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3.2 The Zajc model in the quantum algorithm 

We now interprete the classical probability distribution of ([H|) as a distri¬ 
bution of phase space points associated with the centers of Gaussian wavepack- 
ets according to (0). The one-particle spectrum and two-particle correlator are 
then readily calculated via the emission function (0), 

^ m N _ ( p 2 

Ep d^> ~ (27t) 3 / 2 (Pq + l/2cr 2 ) 3 / 2 6XP V 2P 0 2 + l/a 2 

C(q, K) = 1 + exp ( - q 2 f? 2 m ) , (24) 

= 1 + ^ 2p2 (o>Pg + R 2 J° 2 + 2f? 2 P 2 (l - « 2 )) • (25) 

In contrast to the classical algorithm, the radius parameter is now always posi¬ 
tive, irrespective of the value of s. Even if the classical distribution p^^ s (r, P> *) 
is sharply localized in phase space, its folding with minimum uncertainty wave 
packets leads to a quantum mechanically allowed emission function S^a^p) 
and to a correlator with falls off with increasing q 2 as expected. However, 
the spread of the one-particle momentum spectrum (|2^) receives an addi¬ 
tional contribution l/cr 2 . Choosing a too small increases this term beyond 
the phenomenologically reasonable values, choosing it too large|-|Widens the 
corresponding HBT radius parameters significantly. It was arguedH that a can 
be interpreted as quantum mechanical ‘size’ of the particle, a ~ 1 fm. Given 
the heuristic nature of these arguments and the significant modifications this 
implies for the spectra (|3|) and (p4|), it is however fair to say that presently 
a mainly plays the role of a regulator of unwanted violations of the quantum 
mechanical uncertainty principle while a deeper understanding of its origin in 
the particle production dynamics is still missing. 



3.3 Numerical simulation in the Zajc model 


We have studied the performance of our Bose-Einstein algorithms by generat¬ 
ing with a random number generator sets {(fj, pi, ii)}ie[i,iv m ] of phase space 
points according to the distribution f\.^ s and comparing the num eric al results 
of our algorithms to the analytical expressions of section 3J and |3.2|. 

In Fig. 0(a), the e-dependence of the HBT-radius parameter (|20|) is de¬ 
picted. For fixed bin width e, the approximation of the true HBT radius pa¬ 
rameter by f? 2 iass( e ) seen become better with increasing Pq, in agreement 
with (||^). For the HBT radius obtained from the quantum version of the Zajc 
model and depicted in Fig. 0(b), the situation is both qualitatively and quanti¬ 
tatively different. Now, the HBT radius is always positive, since the Gaussian 











Fig ure 1: Generic properties of the one-dimensional Zajc model, (a) HBT-radius parameter 
(boh of the classical interpretation as a function of the position momentum correlation s. The 
plot shows the HBT radius for different combinations of the model-parameters Rq and Pq, 
and for different sizes of the bin width e needed in numerical implementations, (b) same 
as (a) for the quantum version of the model, eq. (|25|). The dependence of the HBT radius 
on the wavepacket width a is seen clearly, (c) and (d) The two-particle correlator in the 
classical (c) and quantum mechanical (d) version of the model for different sets of model 
parameters. The numerical results are obtained for 50 events of multiplicity 1000, and agree 
very well with the analytical calculations. 


wavepackets take quantum mechanical localization properties automatically 
into account. Also, the s-dependence of the radius is seen to be much smaller, 
since the wavepackets smear out the unphysically strong position-momentum 
correlations present in the classical distribution p^iass- The HBT radius de¬ 
pends not only on the geometrical size Ro , and on the momentum width Pq 
of the source, but also on the wavepacket width a. As seen in Fig. 0(b), this 
wavepacket width affects the HBT radius and its Po-dependence significantly 
for a > Ra¬ 
in Fig. 0(c) and (d) we have presented the two-particle correlation func¬ 
tions corresponding to these HBT-radius parameters for characteristic values 
of the model parameters. The analytical results, obtained by plotting ( ^9|) and 
(UJ) are compared to the results of an application of the event generator algo¬ 
rithms (0) and ([njj) for a Monte Carlo phase space distribution of phase 
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space points. The present plot is obtained for 50 events of multiplicity 1000. 
The differences between analytical and numerical results have their origin in 
statistical fluctuations and become smaller with increasing number of events 
N ev or event multiplicity N m . The good agreement between analytical and 
numerical results in Fig. |l|(c,d) indicates the relatively low statistical require¬ 
ments of our algorithms. The reason is that both algorithms associate to the 
discrete set of generated momenta p, continuous functions of the measured 
momenta pi, P 2 . This smoothens any statistical fluctuation significantly. For 
the classical algorithm, small deviations between numerical and analytical re¬ 
sults are still visible in Fig. [j](c), while the results of the quantum algorithm 
coincide within line width. This can be traced back to the Gaussian prefac¬ 
tor exp (—cr 2 q 2 /2) in eq. (||) which provides an additional smoothening of 
statistical fluctuations not present in the classical algorithm. 
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